##############################
# Media Measurement Matters  #
# Replication Code           #
# Helper Functions           #
##############################

# The following file contains a set of helper functions used to generate plots
# and calculate coefficient values for the paper "Media Measurement Matters:
# Estimating the Persuasive Effects of Partisan Media with Survey and Behavioral
# Data." 

# Score-level overlapping coefficient ----

# Function to normalize scores to range from 0 to 1 (allows for exclusion of NAs)
normalized <- function(x, ...) {(x - min(x, ...)) / (max(x, ...) - min(x, ...))}

# Function to calculate the overlapping coefficient (OC) for a given grouping 
# variable (e.g., partisanship, ideology)
  # var           = revealed preference score for which to calculate overlap (defaults to
  #                 score excluding portal sites)
  # group_var     = group for which to calculate OC (pid = party ID, ideo = ideology,
  #                 med_pref = stated media preferences)
  # values        = unique values for each subgroup 
  #                 - pid: -1 = Democrat, 0 = Independent, 1 = Republican
  #                 - ideo: -1 = Liberal, 0 = Moderate, 1 = Conservative
  #                 - med_pref: strings; MSNBC, Entertainment, Fox
  # out_contrasts = labels for group contrasts (first label = left/right, second
  #                 label = left/middle, third label = right/middle)

score_overlap <- function(var = "score_noportals",
                          group_var = "pid", values = c(-1, 0, 1),
                          out_contrasts = c("Democrat vs. Republican",
                                            "Democrat vs. Independent",
                                            "Republican vs. Independent")){
  
  # Convert srvy to data frame
  srvy <- srvy %>% as.data.frame()
  
  # Normalize score
  srvy$score_norm <- normalized(srvy[, var], na.rm = T)
  
  # Convert scores to radians
  srvy$score_rad <- srvy$score_norm * 2 * pi
  
  # Create datasets by groups
  left <- srvy %>% filter(get(group_var) == values[1] & !is.na(get(var))) %>% 
    select(score_rad) %>% pull()
  mid <- srvy %>% filter(get(group_var) == values[2] & !is.na(get(var))) %>% 
    select(score_rad) %>% pull()
  right <- srvy %>% filter(get(group_var) == values[3] & !is.na(get(var))) %>% 
    select(score_rad) %>% pull()
  
  # Estimate overlap for groups
  leftright_oc <- overlapTrue(densityPlot(left[which(!is.na(left))])$y, 
                              densityPlot(right[which(!is.na(right))])$y)
  leftmid_oc <- overlapTrue(densityPlot(left[which(!is.na(left))])$y, 
                            densityPlot(mid[which(!is.na(mid))])$y)
  rightmid_oc <- overlapTrue(densityPlot(mid[which(!is.na(mid))])$y, 
                             densityPlot(right[which(!is.na(right))])$y)
  
  # Output results
  out <- list(leftright_oc, leftmid_oc, rightmid_oc)
  names(out) <- out_contrasts
  
  return(out)
}

# Score-level distributions ----

# Function to plot distributions of revealed preference scores across groups
  # var        = var to group on (pid = party ID, ideo = ideology, med_pref = stated media pref)
  # var_levels = values corresponding to each subgroup
  # var_labels = string label corresponding to each subgroup value (in order of var_levels)
  # x_var      = revealed preference score to use (e.g., score, score without portals)
  # x_label    = label for x-axis
  # y_label    = label for y-axis
  # df         = data frame to use for analysis (default to the full survey data)
  # x_limits   = range of x-axis (no limits set by default)
  # weights    = whether to incorporate survey weights into plot (T/F)
  # by         = spacing for x-axis breaks (defaults to 0.5)
  # height     = height of y-axis (defaults to 1.5)
  # colors     = colors to use for each subgroup density

score_plots <- function(var, var_labels, var_levels = c(-1, 0, 1), x_var, x_label, 
                        y_label, df = srvy, x_limits = NULL,
                        weights = FALSE, by = 0.5, height = 1.5,
                        colors = c(blue_mit, grey_dark ,red_mit)){
  
  if(weights == FALSE){ # Plot densities by group without including sampling weights
    plot <- df %>% 
      select(all_of(x_var), all_of(var)) %>%
      filter(get(var) %in% var_levels) %>%
      ggplot() + 
      geom_density_ridges(aes(x=get(x_var), 
                              col= factor(get(var), levels = var_levels,
                                          labels = var_labels),
                              fill= factor(get(var), levels = var_levels,
                                           labels = var_labels),
                              y=factor(get(var), levels = var_levels,
                                       labels = var_labels)),
                          lwd = 1.3, 
                          quantile_lines = TRUE, alpha = 0.5,
                          quantile_fun=function(x,...)mean(x), na.rm = T) +
      scale_color_manual(breaks = var_labels, 
                         values = colors,
                         aesthetics = c("color","fill")) + 
      scale_x_continuous(x_label, breaks = seq(-2, 2, by = by),
                         limits = x_limits) + 
      scale_y_discrete(y_label, expand = expansion(add=c(0.2,height))) + 
      theme_minimal() + 
      theme(legend.position = "none",
            axis.text = element_text(size = 14/5 * 3.25),
            axis.title.y = element_text(margin = unit(c(0, 3, 0, 0), "mm")),
            axis.title.x = element_text(margin = unit(c(3, 0, 0, 0), "mm")))
    
  }else{ # Plot densities by group with including sampling weights
    plot_start <- df %>% 
      select(rake_weight, all_of(x_var), all_of(var)) %>%
      filter(get(var) %in% var_levels) %>%
      ggplot() + 
      geom_density_ridges(aes(height =..density..,
                              weight = rake_weight/sum(rake_weight),
                              x=get(x_var), 
                              col= factor(get(var), levels = var_levels,
                                          labels = var_labels),
                              fill= factor(get(var), levels = var_levels,
                                           labels = var_labels),
                              y=factor(get(var), levels = var_levels,
                                       labels = var_labels)),
                          stat="density", lwd = 1.3, 
                          alpha = 0.5, na.rm = T) +
      scale_color_manual(breaks = var_labels, 
                         values = colors,
                         aesthetics = c("color","fill")) + 
      scale_x_continuous(x_label, breaks = seq(-2, 2, by = by),
                         limits = x_limits) + 
      scale_y_discrete(y_label, expand = expansion(add=c(0.2,height))) + 
      theme_minimal() + 
      theme(legend.position = "none",
            axis.text = element_text(size = 14/5 * 3.25),
            axis.title.y = element_text(margin = unit(c(0, 3, 0, 0), "mm")),
            axis.title.x = element_text(margin = unit(c(3, 0, 0, 0), "mm")))
    
    # Extract raw data from plot to add averages
    plot_raw <- ggplot_build(plot_start) %>% purrr::pluck("data", 1)
    
    # Extract densities corresponding to these means
    if(var == "med_pref"){
      density_lines <- df %>% 
        filter(!is.na(get(var))) %>% 
        group_by(get(var)) %>% 
        summarise(x_mean = weighted.mean(get(x_var),
                                         rake_weight, na.rm = T)) %>% 
        mutate(group = c(2, 3, 1)) %>% 
        left_join(plot_raw, on = "group") %>% 
        group_by(group) %>%
        summarise(x_mean = first(x_mean), 
                  density = approx(x, density, first(x_mean))$y, 
                  scale = first(scale), 
                  iscale = first(iscale))
    }else{
      density_lines <- df %>% 
        filter(!is.na(get(var))) %>% 
        group_by(get(var)) %>% 
        summarise(x_mean = weighted.mean(get(x_var),
                                         rake_weight, na.rm = T)) %>% 
        mutate(group = 1:3) %>% 
        left_join(plot_raw, on = "group") %>% 
        group_by(group) %>%
        summarise(x_mean = first(x_mean), 
                  density = approx(x, density, first(x_mean))$y, 
                  scale = first(scale), 
                  iscale = first(iscale))
    }
    
    # Plot weighted results with average ratings appended
    plot <- df %>% 
      select(rake_weight, all_of(x_var), all_of(var)) %>%
      filter(get(var) %in% var_levels) %>%
      ggplot() + 
      geom_segment(data = density_lines, 
                   aes(x = x_mean, y = group, xend = x_mean, 
                       yend = group + density * scale * iscale,
                       col = factor(group, levels = 1:3,
                                    labels = var_labels)), lwd = 1.3) + 
      geom_density_ridges(aes(height =..density..,
                              weight = rake_weight/sum(rake_weight),
                              x=get(x_var), 
                              col= factor(get(var), levels = var_levels,
                                          labels = var_labels),
                              fill= factor(get(var), levels = var_levels,
                                           labels = var_labels),
                              y=factor(get(var), levels = var_levels,
                                       labels = var_labels)),
                          stat="density", lwd = 1.3, 
                          alpha = 0.5, na.rm = T) +
      scale_color_manual(breaks = var_labels, 
                         values = colors,
                         aesthetics = c("color","fill")) + 
      scale_x_continuous(x_label, breaks = seq(-2, 2, by = by),
                         limits = x_limits) + 
      scale_y_discrete(y_label, expand = expansion(add=c(0.2,height))) + 
      theme_minimal() + 
      theme(legend.position = "none",
            axis.text = element_text(size = 14/5 * 3.25),
            axis.title.y = element_text(margin = unit(c(0, 3, 0, 0), "mm")),
            axis.title.x = element_text(margin = unit(c(3, 0, 0, 0), "mm"))) + 
      coord_cartesian(clip = "off")

  }
  
  print(plot)
  return(plot)
}

# Visit-level overlapping coefficient ----

# Function for calculating overlapping coefficient for visit-level plots
  # df            = data frame to use (e.g., with or without portals)
  # bma           = whether to use BMA or Eady at al. scores (default to BMA)
  # group_var     = grouping variable
  # values        = values of grouping variable
  # out_contrasts = contrast names for output

visit_overlap <- function(df = NULL, bma = T,
                          group_var = "pid", values = c(-1, 0, 1),
                          out_contrasts = c("Democrat vs. Republican",
                                            "Democrat vs. Independent",
                                            "Republican vs. Independent")){
  
  # Identify which score to use
  if(bma == T){
    score_use <- "avg_align"
  }else{
    score_use <- "zeta"
  }
  
  df <- df %>% 
    mutate(score_var = get(score_use))
  
  # Normalize score
  df[,"score_norm"] <- normalized(df[, score_use], na.rm = T)
  
  # Convert scores to radians
  df[,"score_rad"] <- df[,"score_norm"]  * 2 * pi
  
  # Create datasets by groups
  left <- df %>% filter(get(group_var) == values[1] & 
                          !is.na(score_var)) %>% as.data.frame()
  mid <- df %>% filter(get(group_var) == values[2] & 
                         !is.na(score_var)) %>% as.data.frame()
  right <- df %>% filter(get(group_var) == values[3] & 
                           !is.na(score_var)) %>% as.data.frame()
  
  # Estimate overlap for groups
  leftright_oc <- overlapTrue(densityPlot(left$score_rad[which(!is.na(left$score_rad))])$y, 
                              densityPlot(right$score_rad[which(!is.na(right$score_rad))])$y)
  leftmid_oc <- overlapTrue(densityPlot(left$score_rad[which(!is.na(left$score_rad))])$y, 
                            densityPlot(mid$score_rad[which(!is.na(mid$score_rad))])$y)
  rightmid_oc <- overlapTrue(densityPlot(mid$score_rad[which(!is.na(mid$score_rad))])$y, 
                             densityPlot(right$score_rad[which(!is.na(right$score_rad))])$y)
  
  # Output results
  out <- list(leftright_oc, leftmid_oc, rightmid_oc)
  names(out) <- out_contrasts
  
  return(out)
}

# Visit-level distributions ----

# Function to plot distributions of visit-level alignment scores across groups
  # var        = var to group on (pid = party ID, ideo = ideology, med_pref = stated media pref)
  # var_levels = values corresponding to each subgroup
  # var_labels = string label corresponding to each subgroup value (in order of var_levels)
  # x_label    = label for x-axis
  # y_label    = label for y-axis
  # df         = data frame to use for analysis
  # x_limits   = range of x-axis (no limits set by default)
  # weights    = whether to incorporate survey weights into plot (T/F)
  # by         = spacing for x-axis breaks (defaults to 0.5)
  # height     = height of y-axis (defaults to 1.5)
  # colors     = colors to use for each subgroup density
  # exemplars  = exemplar sites to show on the plot
  # bma        = whether to use BMA or Eady et al. scores (defaults to BMA)
  # wsj        = whether to include simulated WSJ visits (App. S - defaults to no)

# Function to generate distributions for each subgroup
# Note: weighted plots use this workaround: https://github.com/wilkelab/ggridges/issues/5 
visit_plots <- function(var, var_levels = c(-1, 0, 1), var_labels, x_label, 
                        y_label, df = NULL, x_limits = NULL, weights = FALSE, 
                        by = 0.25, height = 1.5, colors = c(blue_mit, grey_dark, red_mit),
                        exemplars = c("msn.com","cnn.com",
                                      "news.yahoo.com","msnbc.com",
                                      "foxnews.com","nytimes.com"), bma = T, wsj = F){
  
  # Identify which score to use and set height of exemplar labels
  if(bma == T){
    score_use <- "avg_align"
    exemplar_ints <- c(0.6,0.8,0.6,0.8,0.6,0.8)
  }else{
    score_use <- "zeta"
    exemplar_ints <- c(0.45,0.8,0.45,0.8,0.45, 0.8)
  }
  
  # Specify which score to use in the dataframe
  df <- df %>% 
    mutate(score_var = get(score_use))
  
  # Identify location for exemplar sites
  match_visits <- df %>% 
    group_by(domain_recode) %>% 
    summarise(visits = n(), score_var = mean(score_var)) %>%
    drop_na(score_var) %>% 
    as.data.frame()
  
  # Append WSJ in relevant cases
  `%notin%` <- Negate(`%in%`)
  if(wsj == T & "wsj.com" %notin% match_visits$domain_recode){
    match_visits <- match_visits %>% 
      bind_rows(data.frame(domain_recode = "wsj.com", visits = NA, score_var = 0.2754))
  }else{
    match_visits <- match_visits
  }
  
  # Plot results
  if(weights == FALSE){
    plot <- df %>% 
      select(score_var, all_of(var)) %>%
      filter(get(var) %in% var_levels & !is.na(score_var)) %>%
      ggplot() +  
      geom_density_ridges(aes(x=score_var, 
                              col= factor(get(var), levels = var_levels,
                                          labels = var_labels),
                              fill= factor(get(var), levels = var_levels,
                                           labels = var_labels),
                              y=factor(get(var), levels = var_levels,
                                       labels = var_labels)),
                          lwd = 1.3, 
                          quantile_lines = TRUE, alpha = 0.5,
                          quantile_fun=function(x,...)mean(x)) +
      geom_vline(data=subset(match_visits, domain_recode %in% exemplars),
                 aes(xintercept = score_var),lty=2) + 
      geom_text(data=arrange(subset(match_visits, domain_recode %in% exemplars),
                             score_var),
                aes(x=score_var,y= exemplar_ints[1:length(exemplars)],
                    label=paste0(domain_recode,"\n",sprintf("%.2f", round(score_var,2)))),
                hjust=0,nudge_x = 0.01, size = 3.25, lineheight = 0.75) +
      scale_color_manual(breaks = var_labels, 
                         values = colors,
                         aesthetics = c("color","fill")) + 
      scale_x_continuous(x_label, breaks = seq(-4, 4, by = by),
                         limits = x_limits) + 
      scale_y_discrete(y_label, expand = expansion(add=c(0.7,height))) + 
      theme_minimal() + 
      theme(legend.position = "none",
            axis.text = element_text(size = 14/5 * 3.25, color = "black"),
            axis.title.y = element_text(margin = unit(c(0, 3, 0, 0), "mm")),
            axis.title.x = element_text(margin = unit(c(3, 0, 0, 0), "mm")))

  }else{
    plot_start <- df %>% 
      select(score_var, all_of(var), rake_weight) %>%
      filter(get(var) %in% var_levels & !is.na(score_var)) %>%
      ggplot() +  
      geom_density_ridges(aes(height =..density..,
                              weight = rake_weight/sum(rake_weight),
                              x=score_var, 
                              col= factor(get(var), levels = var_levels,
                                          labels = var_labels),
                              fill= factor(get(var), levels = var_levels,
                                           labels = var_labels),
                              y=factor(get(var), levels = var_levels,
                                       labels = var_labels)),
                          stat="density", lwd = 1.3, 
                          alpha = 0.5, na.rm = T) +
      geom_vline(data=subset(match_visits, domain_recode %in% exemplars),
                 aes(xintercept = score_var),lty=2) + 
      geom_text(data=arrange(subset(match_visits, domain_recode %in% exemplars),
                             score_var),
                aes(x=score_var,y=exemplar_ints[1:length(exemplars)],
                    label=paste0(domain_recode,"\n",sprintf("%.2f", round(score_var,2)))),
                hjust=0,nudge_x = 0.01, size = 3.25, lineheight = 0.75) + 	
      scale_color_manual(breaks = var_labels, 
                         values = colors,
                         aesthetics = c("color","fill")) + 
      scale_x_continuous(x_label, breaks = seq(-4, 4, by = by),
                         limits = x_limits) + 
      scale_y_discrete(y_label, expand = expansion(add=c(0.7,height))) + 
      theme_minimal() + 
      theme(legend.position = "none",
            axis.text = element_text(size = 14/5 * 3.25, color = "black"),
            axis.title.y = element_text(margin = unit(c(0, 3, 0, 0), "mm")),
            axis.title.x = element_text(margin = unit(c(3, 0, 0, 0), "mm")))
    
    # Extract raw data from plot to add averages
    plot_raw <- ggplot_build(plot_start) %>% purrr::pluck("data", 1)
    
    # Extract densities corresponding to these means
    if(var == "med_pref"){
      density_lines <- df %>% 
        filter(!is.na(get(var))) %>% 
        group_by(get(var)) %>% 
        summarise(x_mean = weighted.mean(score_var,
                                         rake_weight, na.rm = T)) %>% 
        mutate(group = c(2, 3, 1)) %>% 
        left_join(plot_raw, on = "group") %>% 
        group_by(group) %>%
        summarise(x_mean = first(x_mean), 
                  density = approx(x, density, first(x_mean))$y, 
                  scale = first(scale), 
                  iscale = first(iscale))
    }else{
      density_lines <- df %>% 
        filter(!is.na(get(var))) %>% 
        group_by(get(var)) %>% 
        summarise(x_mean = weighted.mean(score_var,
                                         rake_weight, na.rm = T)) %>% 
        mutate(group = 1:3) %>% 
        left_join(plot_raw, on = "group") %>% 
        group_by(group) %>%
        summarise(x_mean = first(x_mean), 
                  density = approx(x, density, first(x_mean))$y, 
                  scale = first(scale), 
                  iscale = first(iscale))
    }
    
    # Plot weighted results with averages appended
    plot <- df %>% 
      select(score_var, all_of(var), rake_weight) %>%
      filter(get(var) %in% var_levels & !is.na(score_var)) %>%
      ggplot() +  
      geom_segment(data = density_lines, 
                   aes(x = x_mean, y = group, xend = x_mean, 
                       yend = group + density * scale * iscale,
                       col = factor(group, levels = 1:3,
                                    labels = var_labels)), lwd = 1.3) + 
      geom_density_ridges(aes(height =..density..,
                              weight = rake_weight/sum(rake_weight),
                              x=score_var, 
                              col= factor(get(var), levels = var_levels,
                                          labels = var_labels),
                              fill= factor(get(var), levels = var_levels,
                                           labels = var_labels),
                              y=factor(get(var), levels = var_levels,
                                       labels = var_labels)),
                          stat="density", lwd = 1.3, 
                          alpha = 0.5, na.rm = T) +
      geom_vline(data=subset(match_visits, domain_recode %in% exemplars),
                 aes(xintercept = score_var),lty=2) + 
      geom_text(data=arrange(subset(match_visits, domain_recode %in% exemplars),
                             score_var),
                aes(x=score_var,y=exemplar_ints[1:length(exemplars)],
                    label=paste0(domain_recode,"\n",sprintf("%.2f", round(score_var,2)))),
                hjust=0,nudge_x = 0.01, size = 3.25, lineheight = 0.75) + 	
      scale_color_manual(breaks = var_labels, 
                         values = colors,
                         aesthetics = c("color","fill")) + 
      scale_x_continuous(x_label, breaks = seq(-4, 4, by = by),
                         limits = x_limits) + 
      scale_y_discrete(y_label, expand = expansion(add=c(0.7,height))) + 
      theme_minimal() + 
      theme(legend.position = "none",
            axis.text = element_text(size = 14/5 * 3.25, color = "black"),
            axis.title.y = element_text(margin = unit(c(0, 3, 0, 0), "mm")),
            axis.title.x = element_text(margin = unit(c(3, 0, 0, 0), "mm")))
  }
  
  print(plot)
  return(plot)
}

# Persuasion plots ----

# > Generate ranges for binned groups ----

# Function to generate range of scores for each preference group
  # bin_var       = score variable used to segment bins
  # score_version = which revealed preference score to use (e.g., slant, volume)
  # parentheses   = whether to wrap ranges in parentheses 
  # group         = whether to use forced-choice (1) or free-choice group (0)

gen_ranges <- function(bin_var, score_version = "score", parentheses = FALSE,
                       group = 1){
  if(parentheses == TRUE){
    labels <- srvy %>% 
      filter(!is.na(get(bin_var)) & forcedchoice == group) %>% 
      group_by(get(bin_var)) %>% 
      summarise(min = min(get(score_version), na.rm = T),
                max = max(get(score_version), na.rm = T)) %>% 
      mutate(range = paste0("(", round(min, 2), " to ", round(max, 2), ")")) %>% 
      as.data.frame()
  }else{
    labels <- srvy %>% 
      filter(!is.na(get(bin_var)) & forcedchoice == group) %>% 
      group_by(get(bin_var)) %>% 
      summarise(min = min(get(score_version), na.rm = T),
                max = max(get(score_version), na.rm = T)) %>% 
      mutate(range = paste0(round(min, 2), " to ", round(max, 2))) %>% 
      as.data.frame()
  }
  
  out <- labels$range
  
  return(out)
  
}

# > Generate labels for axes ----

# Define function to create y-axis for persuasion plots
  # break_val = maximum label to show on both extremes
  # by_val    = interval for incrementing axis labels

plot_labels <- function(break_val = 0.2, by_val = 0.1){
  # Set the breaks and by for axis
  break_val = break_val
  by = by_val
  
  # Round breaks
  att_vals <- share_vals <- sprintf("%.2f", round(seq(-break_val, break_val, by = by), 2))
  
  # Append labels to attitude axis
  att_vals[1] <- paste0("\n\n\n", att_vals[1], "\n\n(More\nliberal)")
  att_vals[length(att_vals)] <- paste0("(More\nconservative)\n\n", att_vals[length(att_vals)], "\n\n\n")
  
  # Append labels to sharing axis
  share_vals[1] <- paste0("\n\n\n\n", share_vals[1], "\n\n(Less\nwillingness\nto share)")
  share_vals[length(share_vals)] <- paste0("(Greater\nwillingness\nto share)\n\n", 
                                           share_vals[length(share_vals)], "\n\n\n\n")
  
  # Output list of labels
  return(list(att = att_vals, share = share_vals))
}

# > Estimate persuasion overall ----

# Function to calculate unconditional treatment effects
  # weights = whether to incorporate sampling weights into the regression

overall_plot <- function(weights = FALSE){
  
  # Set reference category for treatment indicator
  srvy$article_forced <- factor(srvy$article_forced, ordered = FALSE,
                                levels = c("Entertainment",
                                           "Fox", 
                                           "MSNBC"))
  
  if(isFALSE(weights)){
    # Fox and MSNBC vs. entertainment
    att_mod <- lm_robust(charter_index ~ article_forced, 
                         data = srvy %>% 
                           filter(forcedchoice == 1),
                         se_type = "HC2")
    
    beh_mod <- lm_robust(actions_index ~ article_forced, 
                         data = srvy %>% 
                           filter(forcedchoice == 1),
                         se_type = "HC2")
    
    # Fox vs. MSNBC
    att_mod2 <- lm_robust(charter_index ~ relevel(factor(article_forced),
                                                  ref = "MSNBC"), 
                          data = srvy %>% 
                            filter(forcedchoice == 1),
                          se_type = "HC2")
    
    beh_mod2 <- lm_robust(actions_index ~ relevel(factor(article_forced),
                                                  ref = "MSNBC"), 
                          data = srvy %>% 
                            filter(forcedchoice == 1),
                          se_type = "HC2")
    
  }else{
    # Fox and MSNBC vs. entertainment
    att_mod <- lm_robust(charter_index ~ article_forced, 
                         data = srvy %>% 
                           filter(forcedchoice == 1),
                         weights = rake_weight,
                         se_type = "HC2")
    
    beh_mod <- lm_robust(actions_index ~ article_forced, 
                         data = srvy %>% 
                           filter(forcedchoice == 1),
                         weights = rake_weight,
                         se_type = "HC2")    
    
    # Fox vs. MSNBC
    att_mod2 <- lm_robust(charter_index ~ relevel(factor(article_forced),
                                                  ref = "MSNBC"), 
                          data = srvy %>% 
                            filter(forcedchoice == 1),
                          weights = rake_weight,
                          se_type = "HC2")
    
    beh_mod2 <- lm_robust(actions_index ~ relevel(factor(article_forced),
                                                  ref = "MSNBC"), 
                          data = srvy %>% 
                            filter(forcedchoice == 1),
                          weights = rake_weight,
                          se_type = "HC2")    
  }
  
  # Attitudinal results
  att_out <- as.data.frame(summary(att_mod)$coefficients[2:3, ])
  att_out2 <- as.data.frame(t(summary(att_mod2)$coefficients[3, ]))
  
  att_out <- smartbind(att_out, att_out2)
  
  att_out$outcome <- "Attitudinal Index"
  att_out$id <- c("Fox vs.\nEntertainment", "MSNBC vs.\nEntertainment",
                  "Fox vs.\nMSNBC")
  
  # Behavioral results
  beh_out <- as.data.frame(summary(beh_mod)$coefficients[2:3, ])
  beh_out2 <- as.data.frame(t(summary(beh_mod2)$coefficients[3, ]))
  
  beh_out <- smartbind(beh_out, beh_out2)
  
  beh_out$outcome <- "Sharing Index"
  beh_out$id <- c("Fox vs.\nEntertainment", "MSNBC vs.\nEntertainment",
                  "Fox vs.\nMSNBC")
  
  # Format results
  colnames(att_out) <- colnames(beh_out) <- 
    c("naive", "se", "t", "p", "lower", "upper", "df", "outcome", "id")
  
  out <- rbind.data.frame(att_out, beh_out)
  
  out <- out %>% 
    mutate(min_cilo = lower,
           max_cihi = upper,
           min_cilo90 = naive + qt(0.05, df) * se,
           max_cihi90 = naive + qt(0.95, df) * se) %>% 
    select(id, outcome, naive, min_cilo, max_cihi, min_cilo90, max_cihi90)
  
  out
  
  return(out)
}

# > Estimate persuasion for each subgroup ----

# The function below calculates persuasion estimates when groups are evenly sized
# (e.g., based on terciles or quantiles). group_vsent_plot below calculates these
# same estimates for groups that are unevenly sized (e.g., based on exemplar sites)
  # var     = variable to bin (score, score_noportals, news_consump)
  # nbins   = number of bins (e.g., 3, 4, 5, 6, 50)
  # labels  = labels corresponding to each bin's range (generated using gen_ranges)
  # weights = whether to include weights 

vsent_plot <- function(var = "score_noportals", nbins = 3, labels, weights = FALSE){
  
  # Set the base category for the regression as entertainment
  srvy$article_forced <- factor(srvy$article_forced, ordered = FALSE,
                                levels = c("Entertainment",
                                           "Fox", 
                                           "MSNBC"))
  
  
  # If using weights, bin respondents into weighted vs. unweighted quantiles
  if(weights == FALSE){
    srvy$bin_var <- ntile(srvy[, var], n = nbins)
  }else{
    srvy$bin_var <- weight_bins(score_var = var, nbins = nbins)
  }
  
  # Run the regressions for both outcomes for each preference group
  out <- as.data.frame(do.call("rbind", lapply(1:nbins, function(x){
    # If using weights, add weights into the regression
    if(isFALSE(weights)){
      # Fox and MSNBC vs. entertainment
      att_mod <- lm_robust(charter_index ~ article_forced, 
                           data = srvy %>% 
                             filter(bin_var == x &
                                      forcedchoice == 1),
                           se_type = "HC2")
      
      beh_mod <- lm_robust(actions_index ~ article_forced, 
                           data = srvy %>% 
                             filter(bin_var == x &
                                      forcedchoice == 1),
                           se_type = "HC2")
      
      # Fox vs. MSNBC
      att_mod2 <- lm_robust(charter_index ~ relevel(factor(article_forced),
                                                    ref = "MSNBC"), 
                            data = srvy %>% 
                              filter(bin_var == x &
                                       forcedchoice == 1),
                            se_type = "HC2")
      
      beh_mod2 <- lm_robust(actions_index ~ relevel(factor(article_forced),
                                                    ref = "MSNBC"), 
                            data = srvy %>% 
                              filter(bin_var == x &
                                       forcedchoice == 1),
                            se_type = "HC2")
      
    }else{
      # Fox and MSNBC vs. entertainment
      att_mod <- lm_robust(charter_index ~ article_forced, 
                           data = srvy %>% 
                             filter(bin_var == x &
                                      forcedchoice == 1),
                           weights = rake_weight,
                           se_type = "HC2")
      
      beh_mod <- lm_robust(actions_index ~ article_forced, 
                           data = srvy %>% 
                             filter(bin_var == x &
                                      forcedchoice == 1),
                           weights = rake_weight,
                           se_type = "HC2")    
      
      # Fox vs. MSNBC
      att_mod2 <- lm_robust(charter_index ~ relevel(factor(article_forced),
                                                    ref = "MSNBC"), 
                            data = srvy %>% 
                              filter(bin_var == x &
                                       forcedchoice == 1),
                            weights = rake_weight,
                            se_type = "HC2")
      
      beh_mod2 <- lm_robust(actions_index ~ relevel(factor(article_forced),
                                                    ref = "MSNBC"), 
                            data = srvy %>% 
                              filter(bin_var == x &
                                       forcedchoice == 1),
                            weights = rake_weight,
                            se_type = "HC2")    
    }
    
    # Attitudinal results
    att_out <- as.data.frame(summary(att_mod)$coefficients[2:3, ])
    att_out2 <- as.data.frame(t(summary(att_mod2)$coefficients[3, ]))
    
    att_out <- smartbind(att_out, att_out2)
    
    att_out$outcome <- "Attitudinal Index"
    att_out$id <- c("Fox vs.\nEntertainment", "MSNBC vs.\nEntertainment",
                    "Fox vs.\nMSNBC")
    
    # Behavioral results
    beh_out <- as.data.frame(summary(beh_mod)$coefficients[2:3, ])
    beh_out2 <- as.data.frame(t(summary(beh_mod2)$coefficients[3, ]))
    
    beh_out <- smartbind(beh_out, beh_out2)
    
    beh_out$outcome <- "Sharing Index"
    beh_out$id <- c("Fox vs.\nEntertainment", "MSNBC vs.\nEntertainment",
                    "Fox vs.\nMSNBC")
    
    # Format results
    colnames(att_out) <- colnames(beh_out) <- 
      c("naive", "se", "t", "p", "lower", "upper", "df", "outcome", "id")
    
    out <- rbind.data.frame(att_out, beh_out)
    out$val <- x
    
    out <- out %>% 
      mutate(bin = labels[x],
             min_cilo = lower,
             max_cihi = upper,
             min_cilo90 = naive + qt(0.05, df) * se,
             max_cihi90 = naive + qt(0.95, df) * se) %>% 
      select(bin, id, outcome, naive, p, min_cilo, max_cihi, 
             min_cilo90, max_cihi90, val)
    
    out
  })))
  
  return(out)
}

# Function to calculate persuasion for each (unevenly sized) subgroup
  # var     = variable to bin (score, score_noportals, news_consump)
  # nbins   = number of bins (e.g., 3, 4, 5, 6, 50)
  # labels  = labels corresponding to each bin's range (generated using gen_ranges)
  # weights = whether to include weights 

group_vsent_plot <- function(var = "score_code", nbins = 3, labels, weights = FALSE){
  
  # Set the base category for the regression as entertainment
  srvy$article_forced <- factor(srvy$article_forced, ordered = FALSE,
                                levels = c("Entertainment",
                                           "Fox", 
                                           "MSNBC"))
  
  out <- as.data.frame(do.call("rbind", lapply(1:nbins, function(x){
    if(isFALSE(weights)){
      # Fox and MSNBC vs. entertainment
      att_mod <- lm_robust(charter_index ~ article_forced, 
                           data = srvy %>% 
                             filter(get(var) == x &
                                      forcedchoice == 1),
                           se_type = "HC2")
      
      beh_mod <- lm_robust(actions_index ~ article_forced, 
                           data = srvy %>% 
                             filter(get(var) == x &
                                      forcedchoice == 1),
                           se_type = "HC2")
      
      # Fox vs. MSNBC
      att_mod2 <- lm_robust(charter_index ~ relevel(factor(article_forced),
                                                    ref = "MSNBC"), 
                            data = srvy %>% 
                              filter(get(var) == x &
                                       forcedchoice == 1),
                            se_type = "HC2")
      
      beh_mod2 <- lm_robust(actions_index ~ relevel(factor(article_forced),
                                                    ref = "MSNBC"), 
                            data = srvy %>% 
                              filter(get(var) == x &
                                       forcedchoice == 1),
                            se_type = "HC2")
      
    }else{
      # Fox and MSNBC vs. entertainment
      att_mod <- lm_robust(charter_index ~ article_forced, 
                           data = srvy %>% 
                             filter(get(var) == x &
                                      forcedchoice == 1),
                           weights = rake_weight,
                           se_type = "HC2")
      
      beh_mod <- lm_robust(actions_index ~ article_forced, 
                           data = srvy %>% 
                             filter(get(var) == x &
                                      forcedchoice == 1),
                           weights = rake_weight,
                           se_type = "HC2")    
      
      # Fox vs. MSNBC
      att_mod2 <- lm_robust(charter_index ~ relevel(factor(article_forced),
                                                    ref = "MSNBC"), 
                            data = srvy %>% 
                              filter(get(var) == x &
                                       forcedchoice == 1),
                            weights = rake_weight,
                            se_type = "HC2")
      
      beh_mod2 <- lm_robust(actions_index ~ relevel(factor(article_forced),
                                                    ref = "MSNBC"), 
                            data = srvy %>% 
                              filter(get(var) == x &
                                       forcedchoice == 1),
                            weights = rake_weight,
                            se_type = "HC2")    
    }
    
    # Attitudinal results
    att_out <- as.data.frame(summary(att_mod)$coefficients[2:3, ])
    att_out2 <- as.data.frame(t(summary(att_mod2)$coefficients[3, ]))
    
    att_out <- smartbind(att_out, att_out2)
    
    att_out$outcome <- "Attitudinal Index"
    att_out$id <- c("Fox vs.\nEntertainment", "MSNBC vs.\nEntertainment",
                    "Fox vs.\nMSNBC")
    
    # Behavioral results
    beh_out <- as.data.frame(summary(beh_mod)$coefficients[2:3, ])
    beh_out2 <- as.data.frame(t(summary(beh_mod2)$coefficients[3, ]))
    
    beh_out <- smartbind(beh_out, beh_out2)
    
    beh_out$outcome <- "Sharing Index"
    beh_out$id <- c("Fox vs.\nEntertainment", "MSNBC vs.\nEntertainment",
                    "Fox vs.\nMSNBC")
    
    # Format results
    colnames(att_out) <- colnames(beh_out) <- 
      c("naive", "se", "t", "p", "lower", "upper", "df", "outcome", "id")
    
    out <- rbind.data.frame(att_out, beh_out)
    out$val <- x
    
    out <- out %>% 
      mutate(bin = labels[x],
             min_cilo = lower,
             max_cihi = upper,
             min_cilo90 = naive + qt(0.05, df) * se,
             max_cihi90 = naive + qt(0.95, df) * se) %>% 
      select(bin, id, outcome, naive, p, min_cilo, max_cihi, 
             min_cilo90, max_cihi90, val)
    
    out
  })))
  
  return(out)
}

# > Comparing Fox and MSNBC ----

# Function to estimate treatment effects for Fox vs. MSNBC
  # var      = variable to bin (score, score_noportals, news_consump)
  # pref_val = group on which to filter the data (numeric, typically 1-3)
  # labels   = labels corresponding to each bin's range (generated using gen_ranges)
  # weights  = whether to include weights  

fox_plot <- function(var = "score_code", pref_val, labels, weights = FALSE){
  
  # Set the base category for the regression as entertainment
  srvy$article_forced <- factor(srvy$article_forced, 
                                levels = c("MSNBC", "Fox", "Entertainment"))
  
  if(isFALSE(weights)){ # Without weights
    att_mod <- lm_robust(charter_index ~ article_forced, 
                         data = srvy %>% 
                           filter(get(var) == pref_val &
                                    forcedchoice == 1),
                         se_type = "HC2")
    
    beh_mod <- lm_robust(actions_index ~  article_forced, 
                         data = srvy %>% 
                           filter(get(var) == pref_val &
                                    forcedchoice == 1),
                         se_type = "HC2")
  }else{ # With raked survey weights
    att_mod <- lm_robust(charter_index ~ article_forced, 
                         data = srvy %>% 
                           filter(get(var) == pref_val &
                                    forcedchoice == 1),
                         weights = rake_weight,
                         se_type = "HC2")
    
    beh_mod <- lm_robust(actions_index ~ article_forced, 
                         data = srvy %>% 
                           filter(get(var) == pref_val &
                                    forcedchoice == 1),
                         weights = rake_weight,
                         se_type = "HC2")    
  }
  
  # Attitudinal results
  att_out <- as.data.frame(t(summary(att_mod)$coefficients[2, ]))
  att_out$outcome <- "Attitudinal Index"
  
  # Behavioral results
  beh_out <- as.data.frame(t(summary(beh_mod)$coefficients[2, ]))
  beh_out$outcome <- "Sharing Index"
  
  # Format results
  colnames(att_out) <- colnames(beh_out) <- 
    c("naive", "se", "t", "p", "lower", "upper", "df", "outcome")
  
  out <- rbind.data.frame(att_out, beh_out)
  out$val <- pref_val
  
  out <- out %>% 
    mutate(bin = labels[pref_val],
           min_cilo = lower,
           max_cihi = upper,
           min_cilo90 = naive + qt(0.05, df) * se,
           max_cihi90 = naive + qt(0.95, df) * se) %>% 
    select(bin, outcome, naive, min_cilo, max_cihi,
           min_cilo90, max_cihi90, val)
  
  return(out)
}

# Regression tables ----

# > Comparing partisan media and entertainment ----

# Function to generate stargazer table for individual subgroups
  # var      = variable containing group indicators
  # nbins    = number of categories in group
  # labels   = labels corresponding to each bin's range
  # out_name = label for saving the resulting output
  # weights  = whether to include weights in the regression
  # caption  = caption for table header
  # controls = whether to include control covariates in regression (default to naive ests)

vsent_reg_table <- function(var = "score_code", nbins = 3, labels, out_name,
                            weights = FALSE, caption, controls = FALSE){
  
  # Set order for the experimental conditions
  srvy$article_forced <- factor(srvy$article_forced, ordered = FALSE,
                                levels = c("Entertainment",
                                           "Fox", 
                                           "MSNBC"))
  
  if(isFALSE(controls)){
    att_mods <- lapply(1:nbins, function(x){
      if(weights == FALSE){
        # Fox and MSNBC vs. entertainment
        lm(charter_index ~ article_forced, 
           data = srvy %>% 
             filter(get(var) == x &
                      forcedchoice == 1))
        
      }else{
        # Fox and MSNBC vs. entertainment
        lm(charter_index ~ article_forced, 
           data = srvy %>% 
             filter(get(var) == x &
                      forcedchoice == 1),
           weights = rake_weight)
      }
    })
    
    beh_mods <- lapply(1:nbins, function(x){
      if(weights == FALSE){
        lm(actions_index ~ article_forced, 
           data = srvy %>% 
             filter(get(var) == x &
                      forcedchoice == 1))    
      }else{
        lm(actions_index ~ article_forced, 
           data = srvy %>% 
             filter(get(var) == x &
                      forcedchoice == 1),
           weights = rake_weight)    
      }
    })
    
    out_mods <- c(att_mods, beh_mods)
    stargazer(out_mods, se = starprep(out_mods), p = starprep(out_mods, stat = "p.value"), digits=2,
              column.labels = rep(labels, times = 2), no.space = TRUE,
              column.sep.width = "-7pt", header = FALSE,
              font.size = "footnotesize", style = "apsr",
              table.layout ="=dc-t-s=n",
              dep.var.labels = c("\\textbf{Attitudinal Index}", "\\textbf{Sharing Index}"),
              covariate.labels = c("Fox", "MSNBC"), label = out_name,
              title = caption, notes.append = TRUE,
              notes = c("\\textit{Note}: The reference category is Entertainment. Both dependent variables range from 0 to 1. \\textit{p} values are", 
                        "based on robust standard errors (HC2 variant)."),
              df = F, omit.stat=c("f", "ser"), model.numbers = FALSE,
              out = paste0(table_path, out_name, ".tex"))
    
  }else{
    att_mods <- lapply(1:nbins, function(x){
      if(weights == FALSE){
        # Fox and MSNBC vs. entertainment
        lm(charter_index ~ article_forced + (gender == 1) + race_white + scale(age) + 
             (educ > 3) + scale(int) + 
             factor(pid, levels = c(0, -1, 1)) + factor(ideo, levels = c(0, -1, 1)), 
           data = srvy %>% 
             filter(get(var) == x &
                      forcedchoice == 1))
        
      }else{
        # Fox and MSNBC vs. entertainment
        lm(charter_index ~ article_forced + (gender == 1) + race_white + scale(age) + 
             (educ > 3) + scale(int) + 
             factor(pid, levels = c(0, -1, 1)) + factor(ideo, levels = c(0, -1, 1)), 
           data = srvy %>% 
             filter(get(var) == x &
                      forcedchoice == 1),
           weights = rake_weight)
      }
    })
    
    beh_mods <- lapply(1:nbins, function(x){
      if(weights == FALSE){
        lm(actions_index ~ article_forced + (gender == 1) + race_white + scale(age) + 
             (educ > 3) + scale(int) + 
             factor(pid, levels = c(0, -1, 1)) + factor(ideo, levels = c(0, -1, 1)), 
           data = srvy %>% 
             filter(get(var) == x &
                      forcedchoice == 1))    
      }else{
        lm(actions_index ~ article_forced + (gender == 1) + race_white + scale(age) + 
             (educ > 3) + scale(int) + 
             factor(pid, levels = c(0, -1, 1)) + factor(ideo, levels = c(0, -1, 1)), 
           data = srvy %>% 
             filter(get(var) == x &
                      forcedchoice == 1),
           weights = rake_weight)    
      }
    })
    
    out_mods <- c(att_mods, beh_mods)
    stargazer(out_mods, se = starprep(out_mods), p = starprep(out_mods, stat = "p.value"), digits=2,
              column.labels = rep(labels, times = 2), no.space = TRUE,
              column.sep.width = "-7pt", header = FALSE,
              keep=c("Constant","article_forcedFox", "article_forcedMSNBC"),
              font.size = "footnotesize", style = "apsr",
              table.layout ="=dc-t-s-a=n",
              dep.var.labels = c("\\textbf{Attitudinal Index}", "\\textbf{Sharing Index}"),
              covariate.labels = c("Fox", "MSNBC"), label = out_name,
              title = caption, notes.append = TRUE,
              notes = c("\\textit{Note}: The reference category is Entertainment. Both dependent variables range from 0 to 1. \\textit{p} values are", 
                        "based on robust standard errors (HC2 variant). All models control for age, gender, race, education,",
                        "partisanship, ideology, and attention to national news."),
              add.lines = list(c("Controls?", rep("\\checkmark", nbins * 2))),
              df = F, omit.stat=c("f", "ser"), model.numbers = FALSE,
              out = paste0(table_path, out_name, ".tex"))
  }
}  

# > Comparing Fox and MSNBC ----

# Write function to generate stargazer table for individual score
foxmsnbc_reg_table <- function(var = "score_code", nbins = 3, labels, out_name,
                               weights = FALSE, caption, controls = FALSE){
  # var: variable containing group indicators
  # nbins: number of categories in group
  # labels: labels corresponding to each bin's range
  # weights: whether to include weights in the regression
  
  # Set order for the experimental conditions
  srvy$article_forced <- factor(srvy$article_forced, ordered = FALSE,
                                levels = c("MSNBC", 
                                           "Fox", "Entertainment"))

  if(isFALSE(controls)){
    att_mods <- lapply(1:nbins, function(x){
      if(weights == FALSE){
        # Fox and MSNBC vs. entertainment
        lm(charter_index ~ article_forced, 
           data = srvy %>% 
             filter(get(var) == x & article_forced != "Entertainment" & 
                      forcedchoice == 1))
        
      }else{
        # Fox and MSNBC vs. entertainment
        lm(charter_index ~ article_forced, 
           data = srvy %>% 
             filter(get(var) == x & article_forced != "Entertainment" & 
                      forcedchoice == 1),
           weights = rake_weight)
      }
    })
    
    beh_mods <- lapply(1:nbins, function(x){
      if(weights == FALSE){
        lm(actions_index ~ article_forced, 
           data = srvy %>% 
             filter(get(var) == x & article_forced != "Entertainment" & 
                      forcedchoice == 1))    
      }else{
        lm(actions_index ~ article_forced, 
           data = srvy %>% 
             filter(get(var) == x & article_forced != "Entertainment" & 
                      forcedchoice == 1),
           weights = rake_weight)    
      }
    })
    
    out_mods <- c(att_mods, beh_mods)
    stargazer(out_mods, se = starprep(out_mods), p = starprep(out_mods, stat = "p.value"), digits=2,
              column.labels = rep(labels, times = 2), no.space = TRUE,
              column.sep.width = "-7pt", header = FALSE,
              font.size = "footnotesize", style = "apsr",
              table.layout ="=dc-t-s=n",
              keep=c("Constant","article_forcedFox"),
              dep.var.labels = c("\\textbf{Attitudinal Index}", "\\textbf{Sharing Index}"),
              covariate.labels = c("Fox"), label = out_name,
              title = caption, notes.append = TRUE,
              notes = c("\\textit{Note}: The reference category is MSNBC. Both dependent variables range from 0 to 1. \\textit{p} values are", 
                        "based on robust standard errors (HC2 variant)."),
              df = F, omit.stat=c("f", "ser"), model.numbers = FALSE,
              out = paste0(table_path, out_name, ".tex"))
    
  }else{
    att_mods <- lapply(1:nbins, function(x){
      if(weights == FALSE){
        # Fox and MSNBC vs. entertainment
        lm(charter_index ~ article_forced + (gender == 1) + race_white + scale(age) + 
             (educ > 3) + scale(int) + 
             factor(pid, levels = c(0, -1, 1)) + factor(ideo, levels = c(0, -1, 1)), 
           data = srvy %>% 
             filter(get(var) == x & article_forced != "Entertainment" & 
                      forcedchoice == 1))
        
      }else{
        # Fox and MSNBC vs. entertainment
        lm(charter_index ~ article_forced + (gender == 1) + race_white + scale(age) + 
             (educ > 3) + scale(int) + 
             factor(pid, levels = c(0, -1, 1)) + factor(ideo, levels = c(0, -1, 1)), 
           data = srvy %>% 
             filter(get(var) == x & article_forced != "Entertainment" & 
                      forcedchoice == 1),
           weights = rake_weight)
      }
    })
    
    beh_mods <- lapply(1:nbins, function(x){
      if(weights == FALSE){
        lm(actions_index ~ article_forced + (gender == 1) + race_white + scale(age) + 
             (educ > 3) + scale(int) + 
             factor(pid, levels = c(0, -1, 1)) + factor(ideo, levels = c(0, -1, 1)), 
           data = srvy %>% 
             filter(get(var) == x & article_forced != "Entertainment" & 
                      forcedchoice == 1))    
      }else{
        lm(actions_index ~ article_forced + (gender == 1) + race_white + scale(age) + 
             (educ > 3) + scale(int) + 
             factor(pid, levels = c(0, -1, 1)) + factor(ideo, levels = c(0, -1, 1)), 
           data = srvy %>% 
             filter(get(var) == x & article_forced != "Entertainment" & 
                      forcedchoice == 1),
           weights = rake_weight)    
      }
    })
    
    out_mods <- c(att_mods, beh_mods)
    stargazer(out_mods, se = starprep(out_mods), p = starprep(out_mods, stat = "p.value"), digits=2,
              column.labels = rep(labels, times = 2), no.space = TRUE,
              column.sep.width = "-7pt", header = FALSE,
              keep=c("Constant","article_forcedFox"),
              font.size = "footnotesize", style = "apsr",
              table.layout ="=dc-t-s-a=n",
              dep.var.labels = c("\\textbf{Attitudinal Index}", "\\textbf{Sharing Index}"),
              covariate.labels = c("Fox"), label = out_name,
              title = caption, notes.append = TRUE,
              notes = c("\\textit{Note}: The reference category is MSNBC. Both dependent variables range from 0 to 1. \\textit{p} values are", 
                        "based on robust standard errors (HC2 variant). All models control for age, gender, race, education,",
                        "partisanship, ideology, and attention to national news."),
              add.lines = list(c("Controls?", rep("\\checkmark", nbins * 2))),
              df = F, omit.stat=c("f", "ser"), model.numbers = FALSE,
              out = paste0(table_path, out_name, ".tex"))
  }
}  

# Regression coefficient plots ----

# > Predicting revealed preference scores ----

# Function for generating panels predicting relative slant/volume based on
# party ID, ideology, stated media preferences, and other variables
  # dv        = dependent variable for regression (e.g., score, news_consump)
  # form_comp = component of regression formula to add (e.g., (ideo == 1) + (ideo == -1))
  # label1    = label for first coefficient (e.g., conservative)
  # label2    = label for second coefficient (e.g., liberal)
  # title     = title for the facet (e.g., Party Identification)
  # weights   = whether to incorporate weights into the OLS model
  # scale     = title for the x-axis, based on the DV in use
  # y_low     = lower limit for coefficient values
  # y_high    = upper limit for coefficient values
  # by        = sequence for incrementing axis with coefficient values

panel_plot <- function(dv = "score_noportals", form_comp, label1, label2, title,
                       weights = FALSE, scale = "Average Alignment Score",
                       y_low = -0.3, y_high = 0.3, by = 0.15){
  
  # Set up formula for each regression model
  form <- formula(paste0(dv, " ~ ", form_comp, 
                         " + days_bin + (gender == 1) + race_white + scale(age) + ",
                         "(educ > 3) + scale(polinfo) + scale(int)"))
  
  # Run regression (with or without weights)
  if(weights == FALSE){
    fit <- lm_robust(form, data = srvy, se_type = "HC2")
  }else{
    fit <- lm_robust(form, data = srvy, se_type = "HC2",
                     weights = rake_weight)
  }
  
  # Generate labels for plotting
  attribute_labels <- c(label1, label2, "Consumes News\nDaily", "Male","White", "Age",
                        "Education: College\nDegree or More",
                        "Political\nKnowledge", "Attention to\nNational News")
  
  # Plot results
  coefs_for_plot <- data.frame(summary(fit)$coefficients[-1, c(1,2)])
  
  coefs_for_plot$scale <- scale
  coefs_for_plot$q <- title
  coefs_for_plot$order <- (nrow(summary(fit)$coefficients)-1):1 # for ordering plot correctly
  
  (coef_plot <- ggplot(data = coefs_for_plot, aes(x = Estimate, y = order)) +
      geom_errorbarh(data = coefs_for_plot, aes(y = order, 
                                                xmin = Estimate + qnorm(0.025)*Std..Error, 
                                                xmax = Estimate + qnorm(0.975)*Std..Error), 
                     height=0,colour = 'black') +
      geom_errorbarh(data = coefs_for_plot, aes(y = order, 
                                                xmin = Estimate + qnorm(0.05)*Std..Error, 
                                                xmax = Estimate + qnorm(0.95)*Std..Error), 
                     height=0, lwd=1, colour = 'black') +
      geom_point(colour = 'black', size = 2.25) + 
      xlab("") + 
      scale_x_continuous(limits = c(y_low, y_high),breaks = seq(y_low, y_high, by = by)) +
      scale_y_continuous("",
                         breaks = (nrow(summary(fit)$coefficients)-1):1,
                         labels=attribute_labels) +
      coord_cartesian(xlim=c(y_low, y_high),ylim=c(1, (nrow(summary(fit)$coefficients)-1))) +
      geom_vline(xintercept = 0,size=.5,colour="black",linetype="dotted") +
      theme_bw() +
      facet_grid(. ~ q) +
      theme(axis.text.y = element_text(size=10, angle = 0, hjust = 1, color="black"),
            strip.text.x = element_text(size = 10, face = "bold")) +  
      theme(axis.title.x = element_text(size=12)) )  
  
  return(coef_plot)
}

# > Predicting extremity within preference groups ----

# Function for generating panels predicting extremity based on
# party ID, ideology, stated media preferences, and other variables
  # dv        = dependent variable for regression (default is pref_mod)
  # form_comp = component of regression formula to add (e.g., (ideo == 1) + (ideo == -1))
  # label1    = label for first coefficient (e.g., conservative)
  # label2    = label for second coefficient (e.g., liberal)
  # title     = title for the facet (e.g., Party Identification)

panel_mod <- function(dv = "pref_mod", form_comp, label1, label2 = NULL, title){
  
  # Set up formula for each regression model
  form <- formula(paste0(dv, " ~ ", form_comp, 
                         " + days_bin + (gender == 1) + race_white + scale(age) + ",
                         "(educ > 3) + scale(polinfo) + scale(int)"))
  
  # Run regression
  fit <- lm_robust(form, data = srvy, se_type = "HC2")
  
  # Set up labels for plotting (note that need fewer labels for media preferences vs. PID/ideology)
  if(label1 == "Prefer MSNBC"){
    attribute_labels <- c(label1, "Consumes News\nDaily", "Male","White", "Age",
                          "Education: College\nDegree or More",
                          "Political\nKnowledge", "Attention to\nNational News")
  }else{
    attribute_labels <- c(label1, label2, "Consumes News\nDaily", "Male","White", "Age",
                          "Education: College\nDegree or More",
                          "Political\nKnowledge", "Attention to\nNational News")
  }
  
  # Plot results
  coefs_for_plot <- data.frame(summary(fit)$coefficients[-1, c(1,2)])
  
  coefs_for_plot$scale <- "Extremity of Preferences"
  coefs_for_plot$q <- title
  coefs_for_plot$order <- (nrow(summary(fit)$coefficients)-1):1 # for ordering plot correctly
  
  (coef_plot <- ggplot(data = coefs_for_plot, aes(x = Estimate, y = order)) +
      geom_errorbarh(data = coefs_for_plot, aes(y = order, 
                                                xmin = Estimate + qnorm(0.025)*Std..Error, 
                                                xmax = Estimate + qnorm(0.975)*Std..Error), 
                     height=0,colour = 'black') +
      geom_errorbarh(data = coefs_for_plot, aes(y = order, 
                                                xmin = Estimate + qnorm(0.05)*Std..Error, 
                                                xmax = Estimate + qnorm(0.95)*Std..Error), 
                     height=0, lwd=1, colour = 'black') +
      geom_point(colour = 'black', size = 2.25) + 
      xlab("") + 
      scale_x_continuous(limits = c(-0.2, 0.2),breaks = seq(-0.2, 0.2, by = 0.1)) +
      scale_y_continuous("",
                         breaks = (nrow(summary(fit)$coefficients)-1):1,
                         labels=attribute_labels) +
      coord_cartesian(xlim=c(-0.2, 0.2),ylim=c(1, (nrow(summary(fit)$coefficients)-1))) +
      geom_vline(xintercept = 0,size=.5,colour="black",linetype="dotted") +
      theme_bw() +
      facet_grid(. ~ q) +
      theme(axis.text.y = element_text(size=10, angle = 0, hjust = 1, color="black"),
            strip.text.x = element_text(size = 10, face = "bold")) +  
      theme(axis.title.x = element_text(size=12)) )  
  
  return(coef_plot)
}

# Average ratings ----

# Calculate average ratings by experimental condition within forced-choice group
  # group_var = stated/revealed preference group (e.g., news_consump_bin3, score_code)

avg_df <- function(group_var){
  charter_mean <- srvy %>% 
    filter(forcedchoice == 1 & !is.na(get(group_var))) %>% 
    group_by(article_forced, get(group_var)) %>% 
    summarise(mean = mean(charter_index, na.rm = T), # Group mean
              se = sd(charter_index, na.rm = T)/sqrt(n())) %>% # Group SE
    mutate(upper = mean + qt(0.975, df = n()) * se, 
           lower = mean + qt(0.025, df = n()) * se) %>% # Upper/lower bounds of 95% CI 
    as.data.frame()
  
  action_mean <- srvy %>% 
    filter(forcedchoice == 1 & !is.na(get(group_var))) %>% 
    group_by(article_forced, get(group_var)) %>% 
    summarise(mean = mean(actions_index, na.rm = T), # Group mean
              se = sd(actions_index, na.rm = T)/sqrt(n())) %>% # Group SE
    mutate(upper = mean + qt(0.975, df = n()) * se, 
           lower = mean + qt(0.025, df = n()) * se) %>% # Upper/lower bounds of 95% CI 
    as.data.frame()
  
  colnames(charter_mean) <- colnames(action_mean) <- c("condition", "group", "mean", 
                                                       "se", "upper", "lower")
  
  charter_mean$id <- "Attitudinal Index"
  action_mean$id <- "Sharing Index"
  
  out <- smartbind(charter_mean, action_mean)
  
  return(out)
}

# Plot average ratings by experimental condition
  # group_var = stated/revealed preference group (e.g., news_consump_bin3, score_code)
  # x_labels  = labels for x-axis
  # x_title   = title for x-axis
  # tilt      = angle at which to tilt x-axis labels (default = no tilt)
  # hjust     = justification of x-axis labels (default = centered)

avg_plot <- function(group_var, x_labels, x_title, 
                     tilt = 0, hjust = 0.5){
  df <- avg_df(group_var)
  
  (avg_plot <- ggplot(df, aes(x = factor(group), 
                              y = mean, 
                              fill = factor(condition, 
                                            levels = c("MSNBC",
                                                       "Entertainment",
                                                       "Fox")))) + 
      geom_bar(position = position_dodge(), stat = "identity", width = 0.65,
               color = "black", lwd = 0.25) +
      geom_errorbar(aes(ymin = lower, ymax = upper,
                        group = factor(condition, 
                                       levels = c("MSNBC",
                                                  "Entertainment",
                                                  "Fox"))),
                    width = 0.155, position = position_dodge(0.65), lwd = 1) + 
      geom_point(size = 2.5, position = position_dodge(0.65),
                 show.legend = FALSE) +
      scale_x_discrete(labels = x_labels) + 
      theme_bw() + 
      scale_fill_manual("Condition",values = c(blue_mit, grey_dark, red_mit)) +
      scale_y_continuous(breaks = seq(0, 1, by = 0.25),
                         limits = c(0, 1),
                         expand = c(0, 0),
                         labels = c("More\nliberal\nattitudes", "", "",
                                    "", "More\nconservative\nattitudes"),
                         sec.axis = dup_axis(name="",
                                             breaks=seq(0, 1, by = 0.25),
                                             labels = c("Less\nwillingness\nto share", "", "",
                                                        "", "Greater \nwillingness\nto share"))) +
      coord_cartesian(ylim = c(0, 1)) +
      facet_wrap(~ id, nrow=1) +
      ylab("Average Rating") + xlab(x_title) + 
      theme(plot.title = element_text(hjust = 0.5, size = 16),
            plot.subtitle = element_text(hjust = 0.5, face = "italic", size = 12),
            axis.title.x = element_text(margin = unit(c(3, 0, 0, 0), "mm"), face = "bold",
                                        size = 12, angle = 0, hjust = 0.5),
            axis.title.y = element_blank(),
            legend.title = element_text(hjust = 0.5, face = "bold", size = 12),
            legend.text = element_text(size = 10),
            axis.text.x = element_text(size = 10, color = "black", angle = tilt,
                                       hjust = hjust),
            axis.text.y = element_text(size = 10, color = "black"),
            legend.position = "bottom",
            legend.box = "vertical",
            legend.background = element_blank(),
            legend.box.background = element_rect(colour = "black"),
            strip.text.x = element_text(size = 12)))
  
  print(avg_plot)
  return(avg_plot)
}

